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Abstract 

The work reported in this paper takes place in a general theoretical overview concerning the generation 
of an ensemble time scale. Different algorithms can be des igned to match the particular needs of users and 
the available sets of clocks and time measurements. In all cases however, the statistical treatment of clock 
data requires at least: 

• the definition of an average time scale, 

• the specification of a procedure to optimize the contribution of each clock, 

• the implementation of a filter on each clockfrequency to provide a means ofprediction. 

Here, the comparative study of two time scale algorithms, devised to satisfy different but related require¬ 
ments, is presented. They are ALGOS(BIPM), producing the international reference TAI at the Bureau 
International des Poids et Mesures, and ATI (NIST), generating the real-time time scale ATI at the National 
Institute of Standards and Technology. In each case, the time scale is a weighted average of clock readings, 
but the weight determination and the frequency prediction are different because they are adapted to different 
purposes. 

The possibility of using a mathematical tool, such as the Kalman filter, together with the definition of the 
time scale as a weighted average, is also analysed. Results obtained by simulation are presented. 


INTRODUCTION 


To keep time is to accumulate, without discontinuity, time scale units as close as possible to the SI 
second as defined in 1967 by “the duration of 9192631770 periods of the radiation corresponding to 
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the transition between the two hyperfine levels of the ground state of the caesium atom” M. 

Time laboratories have at their disposal commercial caesium clocks. But physical devices can fail, so 
these laboratories are inevitably led to keep not one, but several clocks which are together treated as 
an ensemble. Clock readings are then combined through an algorithm designed to raise the stability, 
accuracy and reliability of the time scale above the level of performance of any individual clock in the 
ensemble. 

In the design of a time scale algorithm there is no general solution. Rather the fundamental ingredi¬ 
ents should be artfully mixed to match the available time measurements and the needs of the user. 
Some of these ingredients are the definition of an average time scale, the specification of a weighting 
procedure, the determination of a means to predict clock frequencies and the implementation of a 
filter on measurement noise. 

In this paper the key point of the definition of the time scale is highlighted for different algorithms. 
In a first section, we propose the comparative study of two time scale algorithms: ALGOS(BIPM)! 2 !, 
producing the international reference, TAI, at the Bureau International des Poids et Mesures, and 
ATl(NIST)! 3 !, generating the real-time time scale, ATI, at the National Institute of Standards and 
Technology. Though the weight determination and the frequency prediction are different, because 
they are adapted to different purposes, these two algorithms rely on the same definition of the time 
scale. 

In a second section we emphasize the possible use of the Kalman filter for a time scale. This mathe¬ 
matical tool is first briefly presented and then shown as being valuable help in the efficient processing 
of clock data. Here we show with three examples of algorithms based on Kalman filtering! 4, 5, 6] ^ iat 
this technique can be unpowerful for the elaboration of a time scale if an equation of definition is not 
set. Finally, our own view of how to take advantage of the Kalman filter is given together with results 
obtained from simulated clock data. 


Note: In what follows symbols are defined as: 


t 

iu 

MO 


pr- 

Xi(t)=TA - hi(t ): 


N: 

Xij(t) = hj(t) - hi(t) = Xi(t) - Xj(t): 


date of the time scale update, 
clock identification, 
reading of clock Hi at date t 

(this quantity is not directly accessible by experiment), 

weight assigned to clock //<, 

clock II — i time offset from the time scale 

TA under computation 

(this quantity gives user access to the time scale), 
clocks number, 

measurement between clock IIj and clock Hi at date /, 
time interval between two measurement cycles. 


1. COMPARATIVE STUDY OF ALGOS (BIPM) AND ATI (NIST) 

The detailed analysis of the comparison of the two algorithms, ALGOS(BIPM) and ATl(NIST), has 
been published elsewhere^. Here we give only the main features of that study in order to focus on 
the definition of the time scale in these two cases. 
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1-1 ALGOS(BIPM) 


ALGOS(BIPM) produces the international reference, TAI (temps atomique international), at the 
Bureau International des Poids et Mesures (Sevres, France). The requirement here is for extreme 
reliability and long-term stability. To this end, TAI relies on a large number of clocks of different 
types, located in different parts of the world and connected in a network allowing the precise exchange 
of time data. Although measurements are performed at intervals of r=10 days, the definitive update of 
TAI is obtained from an iterative and post-processed procedure which treats, as a whole, two-month 
blocks of data and so ensures long-term stability^ 2 ’ 7 1. One important consequence is that TAI is a 
deferred-time time scale. 

The first step in the establishment of TAI is the computation of a free atomic time scale, EAL (echellc 
atomique libre), obtained as a weighted average of clock readings. TAI is then derived from EAL 
with a frequency steering in order to ensure accuracy. For each date t of the two-month interval of 
computation [to, to + 60 days], EAL is defined as: 

+ '■!(')! 

EAL(t) = -j5-. (1) 

X> 


In this equation, p, is the weight assigned to clock //,, and /i[-(/) is a time correction applied at date 
t to ensure time and frequency continuity of the scale when the weights of clocks or the total number 
of clocks is changed^ 8 !; 


h\(t) = Xi(t 0 ) + B ip ■ (t - t 0 ), ( 2 ) 

here Bi p (t) is the frequency of clock Hj, relative to EAL, predicted for the period to, t. 

From equation (1) and the above notations, we get the system of equations: 

{ N N 

&*.■(*) = EkMW (3) 

*.)(<) = n(d - *,(<) 

The time measurements are chosen to be non-redundant so that the system (3) is deterministic with 
N equations and N unknowns and so is exactly solvable at each date t. The results are the quantities 
Xi(to + nr) with n = 0, 1, 2, 3, 4, 5 and 6 for each clock Hi. Clock Hi frequency B,(f 0 + 60 days) 
for the two-month interval under computation is obtained as the least squares slope of the quantities 

Xi(t 0 + nr). 

The detailed and complete description of the weighting procedure is described elsewhere^. The 
general principle is that the weight assigned to clock Hi is set to be inversely proportional to the 
frequency variance of the clock over six two-month samples. This ensures the long-term stability of 
EAL and allows deweighting for seasonal fluctuation. An upper limit of weight and a system for the 
detection of abnormal behaviour are also in use in ALGOS(BIPM). 
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The frequency prediction is a one step linear prediction^, the supposition being that each clock most 
likely behaves in the present two-month interval as it did in the previous one. This is the optimal 
estimate for random walk frequency modulation, which is the predominant clock noise for two-month 
averaging time. 

1-2. ATI (NIST) 

The ATI time scale, developed at the National Institute of Standards and Technology (Boulder, Co, 
USA) is used for scientific studies. The basic requirement is to provide definitive access to the time 
scale in near real time, with no post-processing or reprocessing. It is an average time scale derived from 
measurements taken from about 10 commercial clocks located on the site. The ATI (NIST) algorithm 
estimates time, adaptative weight and frequency^ for each contributing clock at each measurement 
cycle, at present r = 2 hours. 

The equations for computing the time scale at date t are based on predicted values £,■ and in for the 
time and frequency of each clock. The predicted time difference *,(t) of clock Hi relative to ATI, for 
the date f, involves the time offset Xi(t — r) obtained from the previous computation and the frequency 
Vi(t — T ) estimated at date t — r and predicted for the next r period. This is written as: 


Xi(t) - Xi(t - r) + yi(t - t) ■ r. (4) 

This equation is completely similar to (2) in the description of ALGOS(BIPM) if t — to is set equal to 

T, 

The definition of the time scale itself is written as: 

N 


*.-(«) = 


-*»(«)] 

j =l 

N 


E p , 


(5) 


where x,j(t) is the non-redundant set of time measurements. 
A trivial transformation of (5) leads to: 


- * i ( 0 ] + *«(0 

= — - n - 

E Pj 

so that we get the system of equations: 


YPi*j(t) — ^Pj£j(<) 

Xij(t) = Xi(t) - Xj(t) 


( 6 ) 


(7) 
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This is a system of N equations with N unknowns, equivalent to system (3) for ALGOS(BIPM). 

The analogy in the definition of EAL and ATI time scales is then complete. 

Weights pi appearing in (5), designed to ensure stability, have been determined in the previous compu¬ 
tation at date t — r with an exponential filter over the time deviations between predicted and estimated 
time differences of the last N r periods! 3 !. The time constant N T is usually set at 20 to 30 days. These 
time deviations are also corrected for the bias introduced by the correlation between the clock itself 
and the average time scaled. A detector of abnormal behaviour and an upper limit of weight! 7 ! exist 
also in ATl(NIST). 

The predicted frequency ?/,(/) comes from an exponential weighted average of past and present mean 
frequencies. The time constant of this exponential filter being characteristic of the statistical behaviour 
of each contributing clock. 


Conclusions 

The ALGOS(BIPM) and ATl(NIST) algorithms rely on the same basic definition of the time scale, 
generated as a weighted average of clock readings. They also present other common features: mea¬ 
surements of time differences are treated as having negligible uncertainties and clocks are supposed 
uncorrelated among them. 

The appropriate way to determine clock weights and to predict clock frequencies depends essentially on 
the available measurements (number of clocks, measurement sampling) and on the properties required 
for the resulting time scale (real-time updating or deferred-time post processing). 


2. TIME SCALES BASED ON KALMAN FILTERING 

2-1 KALMAN FILTERING OUTLINE 

The Kalman filter, which is used in many signal processing applications, is a tool well-adapted for 
stochastic estimation and prediction. It is a recursive and linear filter, optimal in the sense of least 
squares estimation! 9, 10 !. 

Its property of recursivity makes of this filter an interesting tool for help in the elaboration of a time 
scale: it allows a definitive treatment at each measurement cycle and provides a means of prediction 
for the next step. The system under estimation, in the case of a time scale, includes clock time offsets 
and clock frequencies. The evolution with time of these quantities can easily be represented by a linear 
model, linearity being a necessary condition for application of the Kalman filter theory. In addition 
measurement noise and correlation among clocks can naturally be inserted in the model whereas it is 
not the case for the two previous algorithms. 

Here is a brief description of how Kalman filtering operates. The intention is to avoid equations which 
can be found elsewhere! 10 ] but rather present basic ideas in the schematic way of Fig. 1 and 2. 

Consider a dynamical system which evolves linearly with time. At date t its state is represented by 
a vector X(t). We wish to estimate this vector using measurements, obtained with a r measurement 
cycle, for times preceding date t. Suppose that the system, state was estimated at date t — r by the 
vector X(t — r/t — r), a quantity which must be read as estimate of X at date t — r knowing 
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all the measurements up to date t — r. This estimate has an error given by a covariance matrix 
T(f — r/t — t) and represented on Fig. 1. 

According to the model of evolution and to the noise of the model given by the covariance matrix 
Q(t), the transition step of the Kalman filter (Fig. 1) allows us to estimate the predicted state of the 
system at date t, knowing all the measurements up to date t - r: the vector X(t/t - t). The error on 
the estimation of this vector is given by the matrix T(t/t - r), which includes the matrix Q(t). This 
error is larger than the error at date t - r, mainly because the model is not perfect (see Fig. 1), 

We now represent a new measure by a vector Z(t ) which is affected by an error given by the matrix 
R(t) and represented on Fig. 1. From the predicted state at date t and this new information, the 
Kalman filter computes a new estimate of the state of the system according to an “update adjustment” 
described in Fig. 2. This new estimate is represented by the vector X(t/t) affected by a covariance 
matrix T(t/t), the trace of which has been minimized. 

The update adjustment of Fig. 2 builds the new estimate from the old one and from the innovation 
weighted by the Kalman gain. The innovation represents the new information contained in the last 
measurement: it is simply the difference between the real measure and a predicted measure, expected 
from the predicted state at t knowing t — r. The Kalman gain K(t) is given by a complex expression 
involving all the errors which affect the system, mainly T, Q and jRl 10 l. Qualitatively, if the mea¬ 
surement Z(t ) is very good, that is, affected by a very small error, the Kalman gain at date t will 
be large so that the new estimate of the system state will largely rely upon the new observation. On 
the contrary, if the measurement is very bad, the adjustment process of the Kalman filter will tend to 
ignore it. 

At last one more detail: the noises which are involved in the Kalman recurrence must be white noises. 


2-2 EXAMPLES OF APPLICATION OF THE KALMAN FILTERING TO THE 
COMPUTATION OF A TIME SCALE 

The first attempt to apply the Kalman filter to the problem of time scales was performed by Tryon 
and Jones in 1982^. Their system is an ensemble of N clocks. The system state X(t) has 2N 
components: the N clock time offsets /i,(/) and the N clock frequencies j/,(t) relative to an ideal time 
scale. The model integrates each clock time and frequency affected by white frequency noise and 
random walk frequency modulation. The measurement vector Z(t) is composed of the (N-l) time 
differences, measured between each clock and the reference clock. The measurement noise is supposed 
to be negligible. The result of the Kalman recursivity is an estimate of how each clock departs from 
an ideal time scale, but it is found that the error of this estimate always increases with time. This 
non-convergence of the covariance matrix T arises from the lack of observability of the system: N 
quantities /i,(t) being estimated from only (N-l) measurements. In this case the Kalman filter is an 
efficient tool for filtering the data noise but, isolated, it does not have the power to build an average 
time scale. 

Another example is the approach developed by Stein^, where a Kalman filter is applied on the time 
measurements z,j(f) themselves, to smooth out the white phase noise. These filtered measurements 
are then used to predict the time offset of a given clock, relative to the ensemble time, in (N-l) different 
ways, each way passing through another clock of the ensemble. The definitive estimate of this time 
offset comes from a weighted average of these different predictions. The weighted average is defined 
by (3) or (5) and computed with a static and robust Kalman filter. 
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The Kalman filter is also used for time scales as a complement to the ATl(NIST) algorithm for 
frequency step detection. This work, proposed by Weiss and Weissert^, utilizes the results aq(i) of 
ATI in order to realize pseudo-measurements of the frequency of each clock relative to the ensemble 
time. The white noise of these pseudo-measurements is filtered and so gives access to the random 
walk component of each frequency and to the variance of this estimation. They are then tested for 
possible step. 


2.3 A NEW PROPOSAL FOR USING THE KALMAN FILTER IN TIME SCALE 
GENERATION 

Here we propose a new approach for using the Kalman fdter in time scale generation and present 
results obtained with simulated clock data. 

We start with the same hypothesis and defining equations as were used for ALGOS(BIPM) and 
ATl(NIST). 

Suppose an ensemble of N clocks, the frequencies of which are uncorrelated. One clock is chosen as 
the reference. Each day N-l time measurements are performed (r = 1 day). Suppose also that the 
white phase noise of the measurements is smoothed out before the main computation of the time scale 
so that it can be treated as negligible. The ensemble time scale is defined by (3) or (5) as: 

f)A[ii(i+.r)T-w,(« + rai 

+ T ) = 1 -—- (8) 

*=i 


with: £i(t + r) = Xi(t ) + #,■(/) • r, similar to (4), 

where is the predicted frequency for the interval [/, t + r]. 

The weight p; and the predicted frequency p,(t), relative to the time scale, of each clock arc chosen 
outside the main computation to ensure the best long-term stability. 

Now we wish to improve the short-term stability of the scale. For this purpose we use N-l Kalman 
filterings, each of them operating on just two clocks, the reference clock Hj and another one chosen 
among the ensemble Hi. This decoupling supposes that the N-l pairs of clocks are uncorrelated, which 
is theoretically not true, as the same reference clock is involved in each pair. However, one can choose 
the least noisy clock as reference and suppose the coupling to be small. Anyway the correlation of the 
N-l pairs can be easily inserted in a Kalman filter operating on all the pairs together. 

For each Kalman filter, the state of the system is composed of a single quantity, the frequency yij(i) 
at date t of clock Hi relative to the reference clock. The model of evolution of the system is written 
as: 


Vi )(^ T T ) — Vij(t) + a ij 1 (9) 

where aij is white noise driving the random walk frequency modulation of the clock. The Q matrix is 
reduced here to the variance of the white noise a,j. 
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Frequency measurements are deduced from time measurements with the equation: 


+ r) = + X<i(t) = yii (t + r) + ft,- ( 10 ) 

and are affected with white frequency modulation /?,j, with variance R. 

The application of the Kalman filter leads to an estimation of the random walk component of the fre¬ 
quency yij(t + r) of clock Hi relative to clock Hj while smoothing out the white frequency modulation. 

Now, the filtered estimate y[j{t + r) of the frequency yij(t + r) can be introduced in the equation of 
definition of the time scale as: 

+ Vi(t) • r] - + Viji 1 + r ) • T ] 

*j(t + T) = —- jf 1 - (11) 

E* 

or 


- *«( 0 ] £[»•■(*)" y'ntt + r ) ■ r ] 

*j(t + r) = — --+ ~-—-— ( 12 ) 

X* E* 

t=i i=i 

The first term of (12) is the time offset of the reference clock Hj with respect to the average time 

scale. The second term is the weighted average of estimations of the frequency of clock Hj relative to 
the time scale, obtained through clock Hi and the filtered frequency of clock Hi relative to clock IIj] 
it is then the frequency of the reference clock Hj relative to the average time scale at date t + r. 

Our proposal consists in filtering the white frequency modulation to estimate the random walk compo¬ 
nent of the frequency of a clock relative to another clock and then introducing this filtered frequency 
in the definition of the average time scale. This approach is thus opposite to that developed by Weiss 
and Weissert^. 

This new procedure has been investigated with simulated clock data: 6 clocks were simulated with 
different levels of white frequency modulation and random walk of frequency for a 300-day period. 
One clock has better short-term and long-term stability than the others: this is chosen as the reference 
clock. The frequency stability for the 5 pairs of clocks is given on Fig. 3. After filtering of the white 
frequency modulation, the short-term frequency stability for each pair of clocks is largely improved, 
as shown on Fig. 4. The efficiency of the filtering is presented on Fig. 5 for a given pair of clocks: the 
white frequency modulation is smoothed out, leading to the extraction of the random walk component 
of the frequency of one of the clocks relative to the other. 

For the computation of the time scale, the weight p; of clock Hi is chosen to be the reciprocal of its 
Allan variance computed over 30 days. The predicted frequency of clock if, relative to the time scale, 
j/,-, is the average of the previous 30-day frequency data. This averaging time is chosen to improve the 
long-term stability of the average time scale. 
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The frequency stability of the resulting time scale, computed either with raw data or after implemen¬ 
tation of the Kalman filtering procedure, is presented on Fig. 6: the average time scale obtained with 
filtered data has a lower level of white frequency modulation and so is more stable using averaging 
time in the range 1-30 days. After a 30-day averaging time, the random walk frequency noise is 
predominant and the two time scales have the same behavior. 

CONCLUSIONS 

The first step of the construction of a time scale is the definition of the ensemble time. For most time 
scale algorithms used in timing centers, the ensemble time is a weighted average of clock readings. The 
determination of the contribution of each clock and the mode of prediction of its frequency relative to 
the time scale are chosen in order to match special user needs and available time measurements. 

The Kalman filter is a tool well adapted to time scale generation once the definition of the ensemble 
time has been given. It helps to smooth out the white phase noise of the time measurements. Its 
use for filtering the white frequency modulation of the clocks themselves is a new approach. In this 
case the time scale is built with the random walk component of the frequencies of the clocks relative 
to a single clock chosen as reference. The short-term stability of the resulting time scale is then 
significantly improved. 
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Figure 1: Recursive procedure for the Kalman filter. 
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Figure 2: Schematic description of the adjustment process for the Kalman 
filter, leading to the updated estimation of the state of the system (The H 
matrix, not specified in the text, is the 'observation matrix' which links 
the vector of the system state, X, to the vector of measurement, Z. If all 
the quantities under estimation are observable, H is the matrix Identity). 
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QUESTIONS AND ANSWERS 


Unidentified Questioner: Was your data tested for periodicity, such as for the apparent ten day 
period in the graph? Was any Chi square test performed? Was the noise tested for whiteness and 
Gaussian behavior? 

Ms. Tavella: Each type of noise was obtained from random white noise. 
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